Definitive Evidence for Order-by-Quantum-Disorder in Er2Ti207 
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Here we establish the systematic existence of a U{1) degeneracy of all symmetry- allowed Hamil- 
tonians quadratic in the spins on the pyrochlore lattice, at the mean- field level. By extracting the 
Hamiltonian of Er2Ti2 07 from inelastic neutron scattering measurements, we then show that the 
[/(l)-degenerate states of Er2Ti207 are its classical ground states, and unambiguously show that 
quantum fluctuations break the degeneracy in a way which is confirmed by experiment. This is the 
first definitive observation of order by disorder in any material. We provide further verifiable con- 
sequences of this phenomenon, and several additional comparisons between theory and experiment. 



Models with frustrated interactions often display an 
"accidental" ground state degeneracy in the classical 
limit. Within mean field theory (MFT), the classical 
degeneracy extends to one of the free energy, even for 
quantum spins. Theoretically, quantum or thermal fluc- 
tuations may lift this degeneracy and thereby select and 
stabilize an ordered state. This phenomenon is called 
"order-by-disorder" (ObD) [T], and has been discussed 
theoretically for more than 3 decades. 

While ObD could therefore be expected to arise fairly 
frequently, it has so far escaped indisputable experi- 
mental detection, to a large extent because of the dif- 
ficulty of distinguishing fluctuation effects from those of 
weak interactions that explicitly break the degeneracy at 
the mean-field level. Hence, to unambiguously identify 
ObD in a material, we need both a detailed knowledge 
of the material's Hamiltonian and a proof that a mean 
field degeneracy exists which is robust to weak perturba- 
tions. We provide both here for the rare earth pyrochlore 
Er2Ti207, and confirm the ObD physics through con- 
frontation of the theoretically-predicted order with ex- 
perimental observations. 

Prior work identified Er2Ti207 as an "XY" antiferro- 
magnet with an ordered ground state [3{8] in zero field. 
ObD was actually already suggested for it [2, 3 , but 
based on an ad-hoc model which led to several signif- 
icant conflicts with experiment, and as such Er2Ti207 
has been regarded as a long-standing puzzle. Our model 
and theory go well beyond this early work and resolve all 
the prior enigmas. Relation to prior work on this mate- 
rial will be returned to at the end of the paper. 

We proceed as follows. First, we prove that, at 
the mean-field level, any symmetry-allowed Hamilto- 
nian for any magnetic material on the pyrochlore lat- 
tice, quadratic in the spins, possesses a U{1) degeneracy, 
which can only be broken by fluctuations or disorder. 
We next extract the parameters of the nearest-neighbor 
model for Er2Ti207 from the fits of linear spin wave the- 
ory with single-crystal high-field inelastic neutron scat- 
tering, show that MFT describes Er2Ti207 well, and that 



the U{1) degeneracy of its model applies to its zero-field 
ordered phase. We then calculate the splitting due to 
quantum fluctuations, and show that the selected state 
is compatible with zero-field measurements. We also pre- 
dict correspondingly a spin-wave gap of ~ 260 mK (and 
other effects) which may be measured in future experi- 
ments. 

General U{1) degeneracy: We project the Hamiltonian 
to that of effective 6* = 1/2 quantum spins describing the 
magnetic doublet of each rare earth ion on the pyrochlore 
lattice. The most general form of H involving two- spin 
interactions is = i • J^^ , where is the /i*^ 
component of the spin on the site i, in the global (x, y, z) 
basis. It is implicit that the symmetries of the pyrochlore 
lattice constrain the relations between the J^^ [9]. The 
mean field (variational) free energy Fmf = Fq^{H — Hq) ^ 
where Hq and Fq are the Hamiltonian and free energy for 
a fiducial system of decoupled spins with applied Zeeman 
fields, is 



+(i + |m,|)ln(i + |m,|)]. 



[(i-|m,|)ln(i-K|) 



(1) 



where P = l/(/c^T), where T is the temperature and 
is Boltzmann's constant, and where = (S^), = 
(iSf) and thus |m^| < 1/2. The entropic part of the 
free energy, i.e. the last term of Eq. ([T]), is obviously 
independent of the orientation of the magnetization m^. 
Now consider the Ansatz 



m^j{a) = pRe e"^" (^kj + ihj^ 



(2) 



where p G [0,1/2], a G [0,27r[, and and hj are the 
local X and y unit vectors, respectively (see Supplemental 
Material), which depend only upon which of the four 
sublattices the site resides. In words, Eq. ([2| describes 
translational invariant states (no unit cell enlargement) 
where all spins make the same angle with their local x- 
axis. (Note that this spin configuration carries no total 



2 



net moment.) This is the manifold of ground states 
identified in Ref. 2 for ErsTisOr. Now, let $ = pe^" = 
+ i^2^ ^i7*^2 ^ Up to an unimportant constant, 
the free energy for the Ansatz Eq. ([T]) as a function of ^ 
reads 

<f[^] =«^' + «*(^*)' + ^l^l', aeC.beR, (3) 

since Eq. ([T]) is quadratic in the spins. Cubic symmetries 
then impose that a = a* = 0, so that depends on 
\^\ only, i.e. solely on |m^|. Indeed, under the three- fold 
rotation along the 111 axis, one finds a ^ a + 27r/3, or 



0, 



(4) 



since Fy^p should remain invariant under the above trans- 
formation. Thus, within MET, the degeneracy is present 
for arbitrary two-spin interactions. Similar arguments 
show that the leading order term splitting the degeneracy 
in the free energy and consistent with cubic symmetry is 



F6 = -c($'' + ($*n, 



(5) 



with some real constant c. Since there is no general ar- 
gument to make c vanish, we conclude that the U{1) de- 
generacy is an artifact of the approximations introduced 
so far. In MET, it is, however, remarkably robust: six 
spin interactions would be required to induce a term of 
the form of Eq. ([5|. In Er2Ti207 (and indeed most other 
rare earth pyrochlores), this is entirely negligible [10 . 
This leaves only fiuctuations - i.e. ObD - to determine 
the splitting coefficient c. 

Local minimum: By expanding about the degener- 
ate states described by Eq. ([2|, we find that for arbi- 
trary (symmetry preserving) exchange parameters, the 
states in Eq. ([2| are extrema of the free energy (see 
Supp. Mat.). Whether or not they are global minima, 
i.e. whether or not they constitute ground states of the 
problem, depends on the parameters J^^-^. We now pro- 
ceed to the extraction of the latter from experiment, and 
lift any potential suspense: for parameters relevant to 
Er2Ti207, these are the lowest-energy states. 

Er2Ti207 Hamiltonian: The effective 6' = 1/2 de- 
scription applies to Er2Ti207 below about 74 K [2) [TT]. 
Nearest-neighbor exchange dominates, for which the 
Hamiltonian is constrained by symmetry to the form [9] 



^ = E [-^-sfsi - J±(s+s- + S-S+) 

(ij) 



./±±biS+s+ + 7*-s-s-] 

■J,± [5l{Qj5+ + QS-)+i^j] 



(6) 



where the sans serif characters Sf denote components of 
the spins in the local pyrochlore bases, where 7 is a 4 x 4 
complex unimodular matrix, and ( = —7* [9 . The linear 
combinations relating the ^(J-^'s (for nearest-neighbor i 
and j) to Jzz^ J±^ Jz± and J±±, the explicit expression 
of 7 and the local bases used in Eq. (|6| are given in the 
Supp. Mat.. 



To determine the four exchange constants and the two 
components of the ^-tensor speciffc to Er2Ti207, we fit 
inelastic neutron scattering data with the structure fac- 
tor obtained from linear spin wave theory in high field 
applied to the Hamiltonian Eq. (|6|. This method was 
described at length in Ref. 9 (esp. in its Appendix 
C). Experiments were carried out on a single crystal of 
Er2Ti207 grown at McMaster University by the fioating 
zone technique Inelastic neutron scattering by the 
time-of-ffight method was performed at the NIST Center 
for Neutron Research using the Disk Chopper Spectrom- 
eter The incident wavelength of 5 A afforded an 
energy resolution of 0.09 meV. Two orientations of the 
crystal were used such that the vertical axes, i.e. the 
crystallographic directions parallel to the applied field, 
were [110] and [111]. Using two field orientations allowed 
an exceptionally comprehensive study of the high-field 
spin- wave spectra. Furthermore, the understanding of 
the zero-field spectra from the ordered state was also 
enhanced by access to the two inequivalent scattering 
planes normal to the field directions. In all color con- 
tour plots herein, the last two panels represent scatter- 
ing within the plane normal to [111]. All others include 
scattering vectors normal to [110]. 

Spin wave spectra arising in the polarized quantum 
paramagnetic state at = 3 T and T = 30 mK were 
fit to the general anisotropic exchange model of Eq. (|6| 
by matching the dispersions in several directions using 
a least squares method. The full 5'(Q,co') was not fit 
to the data, but followed directly from the Hamiltonian 
extracted from the fit to the dispersions. Within the 
linear spin wave approximation and the nearest-neighbor 
model, we find = 2.45 ± 0.23 and g^^y = 5.97 ± 0.08 
f Ref. [141 finds gz = 2.6 and g^y = 6.8), and in 10~^ meV 



J±± 

J zz — 



4.2 ±0.5, 
-2.5 ±1.8, 



J±-- 

Jz± 



6.5 ±0.75, 
= -0.88 ±1.5 



(7) 



Note that these parameters include the nearest-neighbor 
component of the dipolar interactions, and that weaker 
further neighbor components cannot break the U(l) de- 
generacy, as shown above. 

The above parameters Eq. ([t]) place Er2Ti207 in a re- 
gion of the Jzz — J± — Jz± — J±± phase diagram far 
from spin ice. Notably, in sharp contrast to Yb2Ti207 
[9], the interactions J± and J±± involving the local XY 
components of the spins are dominant. Here conven- 
tional magnetic order is expected at low temperature [15] , 
and Curie- Weiss MET is a good starting point. Within 
the latter, we obtain the U{1) degenerate manifold as 
the zero-field ordered states. Other predictions of MET 
compare well with experiment. MET predicts a con- 
tinuous ordering transition at = 2.3 K which im- 
plies a fiuctuation parameter / = /Tc ~ 2.1, given 
the experimental transition temperature Tc = 1.1 K [5]. 
This is much smaller than typical values of / for sys- 
tems with strong quantum fiuctuations (c.f. / = 13 for 
Yb2Ti207 13), and likely largely due to the usual ther- 
mal fiuctuation effects neglected in MET. The zero tem- 
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FIG. 1. The measured S{Q,uj) at T = 30 mK, if = 3 T sliced along several directions. The first five columns show S{Q,uo) 
in the HHL plane, with the field applied along [110], while the last two columns show S{Q,(jj) for the field along [111]. Top 
row: measured S{Q,uj). Bottom row: calculated S{Q,uj), based on an anisotropic exchange model with six free parameters 
(see text) that were extracted by fitting to the measured dispersions. 



perature field-induced transition (for a (110) field) with 
^MF _ 2 74 X, agrees perfectly with the experimental 
value He = 1.7 ± 0.05 T [16]. 

Zero-point fluctuations: Neglecting the tiny six 
spin couplings, only zero-point quantum fluctuations can 
break the degeneracy of a clean crystal at low tempera- 
ture. We show below that they do, though weakly, find 
the preferred states, and quantitatively estimate the en- 
ergy splitting of the degenerate manifold. 

In the spin wave approximation, the energy of the zero- 
point fluctuations per unit cell is given by 

4 r 

^r=^B"zE/ ^kA (8) 

where the sum runs over the four spin wave modes (see 
Ref. ini, and where Vbz is the volume of the Brillouin 
zone. The spectrum 00^. for states as described by Eq. ([2| 
depends on the angle a as illustrated in the Supplemen- 
tal Material, so that eg^ does as well. Performing the 
integration in Eq. ([8| numerically for different values of 
the phase a, we indeed find that zero-point fluctuations 
break the U{1) degeneracy, and that the six equivalent 
values a = nir/S (n = 0, 1, . . . , 5) are the minima of eg^ 
as illustrated in Figure [2] The energy splitting fits well, 
up to a constant, to eg^ = —A/2 cos 6a (c = 32Nu.c.^ in 
Eq. ([5| at T = 0, where Nu.c. is the number of unit cells), 
with A = 3.5 X 10~^meV. The six a = nn/S states are 
equivalent, i.e. related to one another by cubic symme- 
tries, but differ in the absolute orientation of the spins. 
A zero-field cooled sample would be expected to form a 
multi-domain state with an equal volume fraction of each 
state. Indeed, we find that an equal superposition of the 
spectra of all six domains compares well with the exper- 
imental zero field neutron spectrum (see Supp. Mat.). 

Implications: The first prediction of the ObD cal- 
culation is a definite set of six zero-field ground states, 
with a = n7r/3, selected by the positive coefficient A. 
These are exactly the states identified in Ref. 2 , Gen- 
eral symmetry arguments predict either these ip2 states 
or the alternative sequence that would be selected were 



e'^^ (meV) 




FIG. 2. Zero-point fluctuation energy eg^ in the classically 
degenerate manifold parametrized by a. The peak-to-peak 
energy is A 3.5 10~^ meV. 



A < 0, with a = 7r/6 + n7r/3, which are denoted states 
in Ref. 2, The crucial experiment to distinguish the two 
was already noted in this reference: a magnetic field ap- 
plied along (110) to a zero- field cooled sample should 
lead, due to domain alignment, to a sharp increase of the 
(220) Bragg peak intensity for the iIj2 states, but a sharp 
decrease of intensity for the ipi states (see Supplemental 
Material). A sharp increase is consistently observed in 
several experiments Here we make an extensive 
comparison (see Figure [3]) of theory (Supp. Mat.) to 
experimental intensity versus field at Gve Bragg peaks 
including (220), which gives strong evidence for the cor- 
rectness of the ?/^2 ground state and the Hamiltonian pa- 
rameters. The iIj2 state was also found by a sophisticated 
neutron spherical polarimetry study [4 . 

The second consequence of our ObD scenario is the 
existence of a pseudo-Goldstone mode which acquires a 
small gap at low temperature. It is important to empha- 
size that the exchange Hamiltonian in Eq. (|6| has only 
discrete (point group) symmetries, so the appearance of 
a Goldstone-like mode should be surprising! Though no 
surprise seems to be expressed in the literature, the ex- 
istence of such a mode is apparent from multiple reports 
of a large low temperature specific heat |5] ITTHTQ] 
in Er2Ti207. The pseudo-Goldstone mode is also explic- 
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H (in T) 



FIG. 3. Evolution of the Bragg peak intensities with a field 
H II [110]. The experimental data points from Ref. 5 are 
overplotted on the theoretical curves (overall vertical scale 
of experiment was adjusted by hand) obtained when all six 
domains occupy an equal fraction of the volume in zero field. 
The experimental values for (111) and (113) are suppressed by 
instrumental complications, which are partially compensated 
for here by a multiplication factor of 1.3 (see Supp. Mat. for 
more details) . The dashed vertical line shows the critical field 
ifMF ^ -^_74 rp Qb^g^ined within MFT. 



itly visible in our zero field inelastic neutron scattering 
spectra. One can estimate the specific heat by Debye 
theory, Cy = 4:Nu.c. crT^, where Nu.c. is the number of 
unit cells in the system, and 



120 



(9) 



Here a is the usual cubic lattice spacing, and v is the 
geometric mean spin wave velocity (see Supp. Mat.). 
Using the theoretical value for v one obtains dth ^ 
3.6 J • • mol~^. The experimental value from Ref. 5 
(extracted in the Supp. Mat.) is (Texp = 4.6 in the same 
units, comparable with theory. 

Evidently the gap is not visible in current experiments. 
We now estimate it using field theory. Consider the ef- 
fective (Euclidean) action of a system at T = with slow 
space and time variations of the angle a: 

S = I ^^dr ^ {d,^af + ^ {draf - ^ cos 6a 

(10) 

where Vu.c. is the volume of the unit cell, and the param- 
eters T] are obtained from spin wave theory (see Supp. 
Mat.). Expanding the cosine above, we find that the gap 
A to the spin waves is 

A = ^18X/r] = V27A ( J± + J^^/2) ^ 0.02 meV. (11) 



This is below the 0.09 meV resolution of the inelastic 
neutron scattering data reported in Ref. but is cer- 
tainly experimentally accessible. The gap should also 
manifest in a crossover from to activated magnetic 
specific heat for T < A/Zc^ (see Supp. Mat.). A nuclear 
Schottky anomaly below 200 mK [17] makes a direct ob- 
servation challenging, but extrapolation of specific heat 
data from Ref. 5 does suggest a gap of approximately the 
right magnitude (Supp. Mat.). 



From Eq. (10), one may also extract the lengths = 
^/h^^/{lSX), which describe the width of domain walls 
between symmetry-related states. We obtain = 
1.86 a = 18.71 A and ^2 = 2.44 a = 24.55 A for Er2Ti207. 
Confrontation of domain wall theory with experiments 
will be addressed in a future publication. 

Relation to prior theoretical work: Prior theoretical 
work had conjectured the existence of order- by-disorder 
in Er2Ti207, based upon a classical Heisenberg model 
with easy-plane single-ion anisotropy, which exhibits an 
extensive degeneracy very different from the U{1) degen- 
eracy discussed here |2l [3] . This model is microscopically 
inaccurate [7 , and moreover the extensive degeneracy 
obtained within it is not robust. The use of a general 
Hamiltonian, the finding of the robust degeneracy, and 
the extraction of the parameters of Er2Ti2 07 are essen- 
tial ingredients for the new and definitive conclusions we 
draw in this work. 

Discussion: The measurement of the gap via neu- 
trons or thermodynamics is a remaining experimental 
challenge, but higher resolution experiments are needed. 
Neutron scattering data on field-cooled materials which 
are expected to contain single domains, i.e. single a's, 
would allow a wonderful synergy of theory and experi- 
ment and show proof of high control on this interesting 
material. The interesting field evolution of the lineshape 
of the Bragg reflections [5] will be returned to in a fu- 
ture publication. We have achieved a conclusive and de- 
tailed understanding of the magnetism of Er2Ti207, and 
most importantly for the first time shed light on a mate- 
rial where order-by-disorder physics is unambiguously at 
play. 

After completion of this paper, a theoretical preprint 
pQ] appeared, which reaches some of the same conclu- 
sions regarding Er2Ti207. 
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SUPPLEMENTAL MATERIAL 

I. LATTICE AND COUPLINGS 

A. Coordinates 

As usual, the coordinate system of the pyrochlore lat- 
tice is such that there is one "up" tetrahedron centered 



at the origin, with its four corners at 

Ro = |(l,l,l), Ri = I (1,-1,-1), (12) 

R2 = I (-1,1,-1), R3 = I (-1,-1,1), (13) 

where a is the cubic lattice spacing (that of the underly- 
ing FCC lattice). In EraTiaOy, a ^ 10.07 A. 



B. Local bases 

The local cubic bases in which the Hamiltonian Eq. (joj) 
is expressed are the following (aj,bj,ej) bases 



eo = (l,l,l)/\/3 
ei = (1,-1, -l)/v^ 
e2 = (-l,l,-l)/\/3 
^ e3 = (-l,-l,l)/\/3, 



ao = (-2,l,l)/V6 

ai = (-2,-l,-l)/v^ 

a2 = (2,l,-l)/V6 

^ a3 = (2,-l,l)/V6 

(14) 

= X a^, such that spin on sublattice i is = 
S+(a, - ihi)/2 + Sr(a, + ihi)/2 + Sfe,. 
The 4x4 matrix 7 introduced in Eq. ^ is 

f 1 w w'^\ 

1 w'^ w 

w w'^ 1 

\w'^ w 1 J 

where w = e^^*/^ is a third root of unity. 



C. Relations between nearest- neighbor exchange 
constants 
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(15) 



J zz 

J± 
J±± 



-(2Ji- J2 + 2(J3 + 2J4)), 

(2Ji — J2 — J3 — 2J4), 
(Ji+ J2-2J3 + 2J4), 



1 



3V2 



( Ji + J2 + J3 - J4 



(16) 
(17) 
(18) 
(19) 



where Ji , . . , J4 are the matrix elements of the exchange 
matrices Jij between nearest-neighbor sites when the lat- 
ter matrices are expressed in the global (x, y, z) basis. 
Specifically, 



ioi 



— J4 Ji Js 

V-J4 J3 jj 



(20) 



and the other matrices Jij are obtained from Jqi by ap- 
plying the appropriate cubic rotations 
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D. Dipolar interactions 

We may very simply estimate the strength of the 
nearest-neighbor dipolar interactions. We use the fol- 
lowing notation 



47r 



(21) 



where jiB is the Bohr magneton and /io = 47r 10~^ 
the universal magnetic constant. The nearest-neighbor 
distance is r^^j) = where a = 10.07 A 

is the lattice constant |H Mapping Eq. (21) 



to the nearest-neighbor Hamiltonian Eq. ([6|, we get 
Jff = 8.010-3 meV, jf^ = -4.6 10"^ meV, J^*^ = 
3.210-2 meV, and = -3.8 10"^ meV, indicating 

that dipolar and exchange interactions are of the same 
order of magnitude and thus compete. 

We note, again, that further neighbor dipolar inter- 
actions will not break the U{1) degeneracy Eq. ([2|, as 
shown explicitly in the main text. 



E. Relation to prior model 

Previous theoretical work [H [3] proposed the following 
model for Er2Ti207: 

(ij) V i J 

In the D = +oo limit this model corresponds to Jzz = 
Jz± = and J± = J/Q^J±± = J/3 in the language of 
Eq. ^. 



II. ILLUSTRATION OF THE Ts STATES 

Each state in the a = nn/S series is characterized by 
a global basis vector along which each of the four spin 
projection magnitude is largest: 

2 1 1 

ma = ± ^ V^^^^^' ^^^^ 

where (xi,X2,X3) = (x, y,z) is the usual global basis, 
and /i is periodic mod 3. For example, for a = 0, where 
ma{0) = a^, this axis is the x axis, while for a = 7r/3, 
this axis is the z axis: 



fmo(7r/3) = (-l,-l,2)/V6, 
mi(7r/3) = (-l,l,-2)/v^, 
m2(7r/3) = (l,-l,-2)/V6, 

^m3(7r/3) = {1,1,2)/VQ- 



(24) 




FIG. 4. One of the spin states of the U{1) degenerate man- 
ifold. The green arrow shows the direction of the local a^ 
axis. The blue arrow denotes the spin vector for one of the 
Fs states. 





a = a = 7r/6 

FIG. 5. a = and a = 7r/6 spin states 



Each state in the a = 7r/6-hn7r/3 series is characterized 
by a global basis vector along which each of the four spins 
have a zero projection. For example, for a = 7r/6, this 
axis is the y axis: 



fmo(7r/6) = (-l,0,l)/\/2, 
mi(7r/6) = (-1,0, -1)772, 
m2(7r/6) = (1,0,-1)/V2, 

[m3(7r/6) = (l,0,l)/\/2. 



(25) 



III. PROOF OF THE EXISTENCE OF A LOCAL 
EXTREMUM 

Here we prove that the degenerate states described by 
Eq. (§, 



m° {a) = pRe e (^a.j + ih 



(26) 



are local extrema. To do so, we first note that, in gen- 
eral, for a translationally invariant state, the spins can 
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be written 



10 



(27) 



where $ = pe*" = +^^2, ^i,^2,V^j ^ I^, and where 
an ahowed set of is such that the twelve-dimensional 
vectors made of the concatenation of {a^j^, {b^}^ and 
{cl}i are orthogonal to one another for all j = 1, 



i.e. (ag^ ... 



and 



10, 



0, 



= for j, I = 1, .., 10 and 



Now, to prove that the degenerate states constitute 
local extrema, we need only show that the Landau free 
energy around this degenerate manifold does not contain 
terms linear in ^ or tpj. Terms which contain "one" $ or 
i/jj only are readily seen to vanish because they are not 
invariant under time-reversal symmetry. The remaining 
terms (i.e. those that contain $ and one t/jj) have the 
general form 
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(28) 



By applying the lattice symmetry transformations to ^ 
and tjjj^ and requiring that the above term be invariant 
under the latter symmetries, thus imposing constraints 
on fj, we find that fj = / * = for ah j = 1, .., 10. This 
concludes the proof. 



rameter {Jzz-,J± etc.) and ^f-factor independently, keep- 
ing the other fit parameters at their best fit values, and 
assessing visually the range of acceptable fits. These 
ranges were (in meV) 



- 0.09 < < 0.02, 
0.03 < J±± < 0.05, 

while the for the ^-factors 



0.04 < J± < 0.085, 
-0.012 < < 0.08, (29) 



1.8 <gz< 3.2, 5.8 < < 6.3. 



(30) 



Notice that the range of acceptable fits is wider for Jzz 
and Jz± than for J±, J±± and Qxy^ which attests to the 
importance of the XY spin components. To obtain "±" 
type uncertainties as quoted in the text, we somewhat 
arbitrarily took 1/3 of the half- width of the interval ob- 
tained here for each coupling constant. The above pa- 
rameter ranges, and the best fit values, should be viewed 
as a more accurate representation of the acceptable fits. 



V. ZERO-FIELD STRUCTURE FACTOR 

On Figure |8j we show the zero-field inelastic neutron 
scattering data, and the theoretical structure factor ob- 
tained by using the parameters fitted at 3 T, assuming 
that the system is made of the six equally represented 
symmetry-related domains a = nir/S and a = 7r/6+n7r/3 
in rows 2 and 3, respectively. The comparison with the 
a = nn/S series is very good, and we highlight a couple 
features which show that the structure factor obtained 
from the a = nir/S domains compares better than that 
obtained with the a = 7r/6 -\- nn /3 domains. 



VI. BRAGG PEAK INTENSITY 



A. Experimental Bragg peak intensity 



IV. FITS 

Three-dimensional neutron scattering data sets, with 
two dimensions in Q and one in energy transfer, were 
obtained by rotating the single crystal of Er2Ti207 in 
1.5° steps about the vertical axis (corresponding to the 
field direction, either [110] or [111]). Energy vs. Q 
slices through these three-dimensional data sets were 
then made in various directions in the measured Q plane. 

We chose five directions for the H||[110] data set and 
two for the H||[lll] data set, shown in the first five and 
last two columns of Figure [7| respectively. The cut di- 
rections are depicted in Figure |6] We performed least 
squares fits to the extracted spin wave dispersions in the 
= 3 T, T = 30 mK data set, using theoretical val- 
ues obtained from a linear spin wave expansion of the 
Hamiltonian described by Equation ([6| in the main text. 

Uncertainties are non-trivial to estimate in multipa- 
rameter fits. We proceeded by varying each exchange pa- 



Intensities for the five Bragg peak positions shown in 
Fig.[3]of the main text were measured at the NIST Center 
for Neutron Research using Disk Chopper Spectrometer 
using the Disk Chopper Spectrometer (DCS). The data 
was obtained by rotating the crystal about the vertical 
[110] axis in 1.5° steps. The total intensity of each peak 
was summed and the nuclear contribution was subtracted 
using the intensities obtained atT = 2K, i^ = OT. 

The field evolution of the intensity of the (220) Bragg 
peak has been confirmed using the FLEX triple-axis spec- 
trometer at the Helmholtz Zentrum Berlin. The field de- 
pendence of the peak intensity was measured at several 
additional field strengths, and is shown in Fig. |9] This 
data agrees with that obtained from DCS for the (220) 
peak. 

A map of the raw elastic scattering data is shown in 
Fig.|9] The dashed white lines indicate the boundaries for 
the "dark angle", an area of higher neutron absorption 
arising from components of the magnet cryostat rotating 
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[110] direction [-211] direction 

FIG. 6. Representations of scattering planes perpendicular to [110] (i.e. the HHL plane) and [111], showing the FCC Brillouin 
zone boundaries and the corresponding zone centers (labelled in terms of the conventional simple-cubic unit cell). Blue lines 
indicate the directions of the seven different cuts shown in Figures [l] [7| and [S] 




(HHH) (OOL) (22L) (HH2) (-H+1 , -H+1 , H+2) (-2H, H+1 , H-1 ) (H-1 , 2, -H-1 ) 

FIG. 7. The dispersion curves, which were fit to the if = 3 T data set, are shown in white, overplotted on the data. 




(HHH) (OOL) (22L) (HH2) (-H+1 , -H+1 , H+2) (-2H, H+1 , H-1 ) (H-1 , 2, -H-1 ) 

FIG. 8. Comparison of the if = T, T = 30 mK data (top row) to the calculation involving the a — ni^ j?) (middle row) and 
a — n/Q + mi/^ ground states (bottom row). In each case, the domains are assumed to occupy an equal fraction of the volume. 
Note that the data in HHH and OOL panels has been corrected for the low Q self- absorption (see Fig. [9| in order to make 
easier visual comparisons to the calculations. Circles indicate regions where the two sets of calculations can be distinguished 
and compared to the data. In particular, for the HHH cut, the a = 7r/6 + n7r/3 series shows relatively strong intensity at H=l, 
when the data displays weak intensity at this point. The maximum observed around H=2 in the a — mi/'^ series matches 
the data much better. In the OOL cut, the situation is similar: the intensity maximum at L=2 in the a = 7r/6 + n7r/3 series 
disagrees with the data, which exhibits a maximum at L=4, like the spin waves of the a — n7r/3 series. In the HH2 cut, the 
maximum of the intensity at H=0 in the a — 7r/6 + n7r/3 series is incompatible with the data, whose intensity maximum at 
H?^l agrees better with the theoretical structure factor of the a — n7r/3 series. We are unable to identify obvious differing 
features on the other cuts. 
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FIG. 9. Raw elastic scattering data at if = 3 T, showing 
the locations of the Bragg peaks. Note that (113) and (111) 
fall within the "dark angle" as described in Section |VI| and 
as denoted by dashed white lines here. The self- absorption 
affecting the low Q region (inside the solid white line) is re- 
sponsible for the unphysically low intensities near Q = in 
the (HHH) and (OOL) slices shown through this paper. 



into the incident beam. This dark angle covers both the 
(111) and (113) Bragg positions, leading to artificially 
low integrated intensities for these peaks. One can cor- 
rect for this by measuring the reduction in incoherent 
elastic scattering in this area, which we have done using 
the high field data set in order to avoid any diffuse mag- 
netic scattering that is present in low fields, and find a 
factor of 1.3 ±0.04 reduction in intensity. Thus, in Fig.jS] 
of the main text, the intensities of the (111) and (113) 
peaks have been multiplied by a factor of 1.3. It should 
also be noted that self-absorption effects are also present 
at low HQ II, and begin to be important at the wave vec- 
tors indicated by the solid white line in Fig.[9j The (111) 
peak thus suffers from an additional decrease in inten- 
sity, which we have not corrected for due to the difficulty 
in accurately doing so. It should also be noted that the 
same self-absorption is responsible for the suppressed in- 
tensity near Q = 0, which is evident in all (HHH) and 
(OOL) slices. 



B. Reminder of the linear spin wave theory of 
Ref. 9 



In this paper, we have made extensive use of spin wave 
theory. It is described in detail in Appendix C of Ref. 9, 
so that we use the same notations and only give here 
those notations and definitions needed to understand the 
calculations performed in this paper. 

The unit vector u^, a = 0, .., 3 is defined to be the di- 
rection of Sa which minimizes the classical energy, and 
and Wa are such that (u^, v^, w^) forms an orthonormal 
basis. We further define the Holstein-Primakoff trans- 
verse bosonic operators Xa = x\ and i/a = yl at each site 



such that [xa, Ha] = ^ and 

Sa - Ua = S - Ua, Sa ' = ^a, Sa ' = Ua, 

(31) 

where s = 1/2, and Ua = {x^^ + — l) /2 measures the 
magnetic moment due to quantum fluctuations repre- 
sented by Xa and ya- 



C. Elastic structure factor at Bragg peaks 

The inelastic structure factor is proportional to 
X(k,a;)= (32) 
^ [5,, - (k)^(k),] ^(M,^(-k,-c^)Mak,c.)), 

a, 6 

where = ^^Qa^S^ is the magnetic moment opera- 
tor. By definition, 

M^^(k,a;) = J dre'^^e'^^M^{k)e-'^^, (33) 

so, defining the above expectation value, at zero temper- 
ature, 



(34) 



we can rewrite, using the usual spectral decomposition, 

MT^iKuj) (35) 
= ^S{co + en-eo) (0|M„^(-k)|n)(n|Mak)|0). 

n 

Here the sum on n runs on the eigenstates of H. We 
define the amplitude 



Ak = (0| ^Ma(k) O) = • (0|Sa(k)|0) (36) 

a=0 a=0 

= Ef J-(OK|0)) ga-u.e^'^i^", (37) 



a=0 



where Ua is the direction of Sa which minimizes the clas- 



sical energy, Ua is as defined in Section [VI B[ and because 
(xa) = iVa) = 0- We obtain the elastic structure factor 
as the zero frequency limit of the inelastic one (or more 
properly, integrating the latter over a narrow interval of 
frequency near zero): 

X(k,a; = 0)(x (38) 
^[V-(k)^(k).] (0|^M,^(-k)|0)(0|^M,^(k)|0) 



= A_k • Ak - (k • A_k) (k • Ak) 
= k X ^k X Ak) , 
since A_k = A^. 



(39) 
(40) 
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It now suffices to obtain Ai^. However, we must ac- 
count for the domain structure. In zero field, the six 
a = utt/S states are the ground states, and we assume 
that they are equally present in the system. Then, the 
average intensity of a Bragg peak at k = Q is 



Jq(f = o) = ^X]|Aq(W3)I^ 



(41) 



n=0 



as shown on Figure 10 In an applied [110] field, there are 
only two equilibrium domains (or one above the critical 
field), and both domains have equal intensity for all the 
Q studied, so no averaging is necessary. 

For certain Bragg peaks, in particular Q = ^ (2, 2, 0), 
which we will also denote Q = 220 in reciprocal lattice 
unit vectors, the alignment of domains leads to a jump 
in intensity. For this Q = 220 peak, the zero field inten- 
sities of the six domains take two large and four small 
values (Fig. 10), while an infinitesimal field selects the 
two domains with large intensities. This leads, ideally, 
to a jump in intensity of a factor of two in passing from 
the six domain to two domain state. 

We note that were we to choose the other sign of A, 
such that the i/ji states were ground states, we would ob- 
tain four large and two small intensity contributions in 
zero field, and a consequent decrease in intensity with 
applied field. This contrasting behavior provides a con- 
clusive proof that the selected states are the ?/^2 ones 
{a = utt/S) and not the t/ji ones (a = 7r/6 + n7r/3). 

The domain averaging is taken in account in Fig. 3 of 
the main text, which gives a similar (but opposite) effect 
for the Q = 113 peak, also seen experimentally. 

|A(e=220)|2 




FIG. 10. Evolution of the Q = 220 Bragg peak intensity 
in zero field, as a function of a. In blue the amplitude of 
the zero-temperature classical state, i.e. when rii = 0, with 
analytic formula | Ao (Q = 220, a) |^ = ^9xy sin^ [a + 7r/6] 
is drawn. In purple, the "fit" | A (Q = 220, a) |^ = 
sin^ [a + 7v/6] /\A (Q = 220, a = 7r/3) |^ The points on 
the latter curve are the points with a = n7r/6, n = 0,..,ll. 



VII. CALCULATION OF THE GAP AND 
RELATED QUANTITIES 

To calculate the gap, we consider H = and a = 
0. Then, classically and at zero temperature, using the 



notations of Ref. 9 , which are also given in Section [VI B[ 
2Sa = Ua = a^, and we can choose = and = 
to form the orthonormal basis (u^, v^, w^) = (a^, 6^,6^). 
With the Ansatz Eq. ([2| in mind, we define 



= (i 1 1 1) , *r = (i -1 o) , (42) 

M/^ = (o 1 -l) , = (l 1 -1 -l) (43) 



and parameterize 



Xo Xi X2 X3 ^ 



a*g^ + ^Xi*r, (44) 



and use, as in Ref. [H = (^y^ y2 y^^. The action 
of the linear spin wave theory developed in Ref. |9] is 



2/3 



S ■ 
where 




^2C^ + i{i0Jn)h 



2Ck - iiiujn)h 
2Bk 



(45) 



(46) 



with Ak, 5k and Ck as defined in Ref. [H If we integrate 
out Y, we are left with 



2(1 



(47) 



where 



with 



iVk.a;„ = [Tk + i(iu;„)li - (iu;„)'Wk] , (48) 



Tk = 2 (Ak — Ok • -L-k ' '-'k 



Vk 



Ck • -Bk ^ • 

(-^k ^ • C'k - Ck ■ B~^) , 



(49) 
(50) 

(51) 



2 ■ 

We want to expand N\^^ojn second order in k and ujn 
about k = and uJn = (we are looking at the pseudo- 
Goldstone mode), so we expand Tk, Vk and to second, 
first and zeroth order, respectively. Then, defining L^"^^ 
as the n}^ order term of matrix L, N\^^^^ is 



A^k,c 



because T^^^ = V^'> = since the k dependence comes 
solely from cos (k • (r^ — r^)) terms, whose expansion in- 
volves even powers of only. Importantly, we fur- 
ther find ■ V(o) • = (and in fact V^^^ = 0), 



.(1) 
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^'J' • T(o) • % = and *J • T^°^ -^j^O for j = 1, 2, 3. 
So, if we complete the squares of the Xj terms, only the 
following terms involve a and Xj- 



2aL 



-a 



(53) 



with ai^^ 
and b^. 



AO) 

~ ^3 



ll.AJ'^Kk. - w'^'\iujn)\ (with 

^ ^) ^^^^ ^^^^^ 
that arise from integrating out Xj ^^'^ \ed.^i of order 

four in k and ujn- This result is reasonable since only 
a describes the continuous degeneracy and is thus ex- 
pected, alone, to give rise to the Goldstone mode. 
We find, for the ground state a = 0, 



S" = 



where 



2/3 



(54) 



X [K.xkl + i^yz{kl + kl) - r]{iujnf 



(2J± - J±±) 



- (4J± + J±±) 
4 1 



3 2 J± + J, 



(55) 
(56) 
(57) 



SO that the action of the full spin-wave and zero-point 
fiuctuation problem is 



S' 



2p 



^ — k,— ^k,a 



(58) 



X [K.^kl + K.y^{kl + kl) - r]{iujnf + 18A] , 



where we expanded the cosine, — | cos 6a — | + 9Aa^ 
and took into account the overall 1/2 prefactor. The gap 
A is then 




27A J± 



Jz 



: 0.0222 meV, (59) 



since A = 3.51 x 10""^ meV. In Kelvins, this is A = 258 
mK. Now, if we define 



we get 

Vr ^ 0.0413 meV and 



: 0.0541 meV 



(60) 



(61) 



where a is the lattice constant, and Vi = Vi/a^s has the 
dimension of an energy. Now, plugging in a = 10.04 A, 
we get 



0.416 meV.A and 



yyz 



0.545 meV.A. 

(62) 



Note that the specific form of the anisotropy, i.e. Vx 7^ 
Vy = Vz is due to the choice a = 0. The other combina- 
tions are found for other values of a = n7r/3. For all the 
latter the velocity which we denote vi appears once, while 
V2 appears twice. Here vi = and V2 = Vy = = Vy^z- 
The read-off slopes of the Goldstone modes of the spin 
wave theory are 



0.0412 meV 



and 



0.0541 meV, 



(63) 

i.e. a basically exact match. We can also define two 
length scales. 



18A' 



(64) 



and we get 



1.86 o w 18.71 A and = 2.44 a « 24.55 A, 

(65) 

where a is, again, the lattice spacing. Those length scales 
physically represent the lengths over which the system 
sees no degeneracy breaking (cf. A = ^ ^ oo), and 
are the typical extent of domain walls in the system, if 
any. 



VIII. SPECIFIC HEAT 




FIG. 11. Cv/T'^ versus T. In blue, the theoretical plot 
Eq. ( 71 ) obtained for the parameters of our fits in zero 
field, with A = 0.02 meV. Above 0.05 K the Unear behav- 



ior of Cv /T'^ is clear, with a slope ath 



Ma 



3.62 J.K~^.mol~^. The purple points are the experimental 
data points reported in Ref. |5] (note that in the latter ref- 
erence, the vertical axes is Cp/R per mole of formula unit 
Er2Ti2 07), and the solid purple line is the best linear fit to 
the six data points with lowest T. The equation of this line is 
-0.2+4.6 TJ.K-^mor^ (aexp = 4.6 J.K-^.moP M. The thin 



dashed line's equation is 



167r4/15 



Qh2 



where r is the ratio of ath to the slope of the the blue line. 
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A. Calculation of the specific heat 



For ease of notation, we take Vi = and V2 



V. 



y 



= '^yz^ which corresponds to the a = state (as de- 
fined and mentioned in Section VII). Of course none of 



the conclusions drawn here depend on this specific choice. 
The dispersion relation of the low-energy spin- wave mode 
is {h = 1) 



where X = /3k and 5 = /3A are dimensionless. 
Cv /{^Nu.c.) is plotted in blue on Figure 



11 



Estimate of the coefficient of the T term 



= ek = ^Jvlkl + v^, {kl + kl) + A2, (66) 

hence (the counting goes: there are as many such modes 
as there are unit cells), following Debye, the energy of 
the system is 



1 



(fke\^- 



1 



BZ 



./3ek 



-y, (67) 



where Vbz = is the volume of the Brillouin zone 

and A^BZ = Nu.c. is the number of points in the Brillouin 
zone (equal to the number of unit cells) . We now change 
the integration variables (rescale) 



kx 



ky^z'^yz 1 



k = ^^^x^, (68) 



The theoretical (blue) curve on Figure 11 is made of 



two parts. Below T > 0.05 K the behavior is that of an 
activated Cy, while the straight line above T > 0.05 K 
clearly pertains to the behavior. Setting A = 0, i.e. 
(5 = in Eq. (71), we extract the slope of this line: 



where 



k% TT^ 



^Nu.c.(tT\ (72) 



3.7510-^ meV.K-^ (73) 



(note k^ has the units of an energy), to get 



E 



1 Nu.c 

Vbz J"bzx^3. 



47r Nu.^ 
Vbz Jo 



dkk^ 



(69) 
(70) 



{vivl 



.1/3 



where = \/k^ + A^, v = {vxVyVz)^^^ 
the geometric mean of the velocities, where the integra- 
tion runs to infinity because we have assumed ^ ksT^ 
where A is the ultraviolet cut-off, A ~ 1/a. The specific 
heat C 



dk ^ 




(71) 



which means, per Er spin. 



Er 



(74) 



with the above value of <j, where Ney is the number of Er 
spins or, in units more commonly used in the literature. 



fiA=0 

— ^ ^ = aT^ with a ^ 3.62 J.K-^mor^, 
moles oi Er 

(75) 

since a = (jJVa (Joules per 1 meV), where A/a is the 
Avogadro constant and (Joules per 1 meV) ^ 1.6010"^^ 
J/meV. Note that a is denoted a in the main text for 
cosmetic reasons. 



